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We study the weak gravitational lensing effects caused by a stochastic distribution of dark matter 
halos. We develop a simple approach to calculate the magnification probability distribution function 
which allows us to easily compute the magnitude bias and dispersion for an arbitrary data sample 
and a given universe model. As an application we consider the effects of single-mass large-scale 
cosmic inhomogeneities (M ^ 10^^h~^ Mq) to the SNe magnitude-redshift relation, and conclude 
that such structures could bias the PDF enough to affect the extraction of cosmological parameters 
from the limited size of present-day SNe data samples. We also release turboGL [1^, a simple and 
very fast (< Is) Mathematica code based on the method here presented. 
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I. INTRODUCTION 

Large-scale inhomogeneities are known to affect the 
light coming from very distant objects. It is important 
to understand these effects well if one is to use cosmologi- 
cal observations to accurately map the expansion history 
and determine the precise composition of the universe. In 
particular the evidence for dark energy in the current cos- 
molgical concordance model is heavily based on the anal- 
ysis of the apparent magnitudes of distant type la super- 
novae (SNe) [2113]. One way inhomogeneities can affect 
the observed SNe magnitude-redshift relation is through 
gravitational lensing. How large these effects are depends 
strongly on the assumptions regarding size and density 
contrast of the structures through which light passes on 
its way from source to observer. The main purpose of this 
paper is to develop a simple tool to compute weak lens- 
ing effects on the SNe magnitude-redshift relation due to 
statistical distributions of large-scale inhomogeneities. 

The effect of inhomogeneities on cosmolgical observ- 
ables have been studied earlier by many authors and in 
many different contexts. The effect of matter clumping 
into isolated halos, or isolated cores, was already con- 
sidered by Kantowski in 1969 [4 , and more recent anal- 
yses of gravitational lensing by statistically distributed 
inhomogeneities have been carried out, for example, in 
Refs. [ail[7lliia[IOl[IIl|lllI3[ll[Ti[T^ Lately 
many authors have studied exactly solvable models for 
large-scale inhomogeneities, such as swiss-cheese, onion 
and meatball models [i[ia[l8l[l9l|20l[ail22l[23lia 
The local Hubble bubble scenario has also been quite 
succesfull in explaining many of the cosmological obser- 
vations m [13 [H m |30l [SH [321 33 , 34, 35 , obvi- 
ously with the price of giving up the cosmological princi- 
ple. Yet another aspect of inhomogeneities is to cause 
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possibly strong backreaction effects on the dynamical 
Einstein equations governing the evolution of the met- 
ric [Ml [33 [Ml [3i [401 [111 [11 [43 [11 [43^ 

The approach in this paper is closest in spirit to that 
of Refs. [6l [12]. Here the cosmological principle is re- 
spected and we will assume that the inhomogeneities can 
be introduced as a perturbation on a well defined global 
background solution. We will also neglect all redshift ef- 
fecs in voids (see Refs. 1171 [H HSl EI]) and concentrate 
purely on cumulative weak lensing [52] . The central phys- 
ical concept for the present work is the observation, made 
already by Zel'dovich in 1964 [53|, that light travelling in 
empty parts of a clumpy but globally homogenous uni- 
verse becomes demagnified. Light going through mass 
concentrations is magnified on the other hand and the 
null result of the Friedmann-Lemaitre- Robertson- Walker 
(FLRW) background solution only arises through aver- 
aging over a large number of photon geodesies. As a re- 
sult, the magnification probability distribution function 
(PDF) for a single observation is skewed favouring mild 
demagnifications, with a long compensating tail of pos- 
itive magnifications. Such skewness is a concern for the 
interpretation of the cosmological data, because the de- 
magnification effect could be misinterpretated as a rela- 
tive dimming of standard candles in small data sets. The 
degree of skewness of the PDF obviously depends on the 
size and the spatial distribution of the matter concentra- 
tions. 

The central object to compute then is the magnifica- 
tion probability distribution function for a random pho- 
ton geodesic in a given inhomogenous universe model. 
Refs. ini [12] evaluated the PDF numerically perform- 
ing large simulations in a model universe and comput- 
ing the magnification factors with ray-tracing techniques. 
Here we will introduce a much simpler way to compute 
the PDF and give an explicit expression for it as a sum 
over probability weighted configurations of inhomogeni- 
ties. Moreover, we will develop a simple analytic approx- 
imation that can reproduce the mode and the dispersion 
of the numerical PDF. In addition to the statistical bias 
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due to skewness, there may be systematic biases on the 
observed PDF, such as extinction, foreground hght con- 
tamination, strong lensing and outher corrections. These 
biases can be problematic because they, unhke the statis- 
tical magnification bias, persist even in large data sets. 
We will not try to estimate the size the systematic bi- 
ases here, but we will show how their effects can easily 
be included in the analysis. 

We test our method by reproducing the key results 
of Refs. [6", ^V2} for a universe in which all matter is ho- 
mogenously distributed in halos with M ~ 10^^ /^"^M© 
and our results are found to be in good agreement. We 
also consider much larger halos of size M ~ 10^^ /^"^M©, 
whose existence is suggested by the large voids and fil- 
amentary structures seen both in the large-scale simula- 
tions and in the galaxy redshift surveys [54l HH [56] . From 
the weak lensing point of view such halos can be consid- 
ered as localized lenses, irrespectively of whether they 
are gravitationally bound or not. We again produce the 
simulated PDFs for a universe with such large structures, 
both in the ACDM and in the Einstein-de Sitter (EdS) 
background model. We also produce the distributions 
for binned sets of observations and compute bias and 
dispersion for these effective PDFs. Quite interestingly, 
supercluster size halos with M ~ lO^^/i^^M© turn out to 
represent roughly the borderline for when the skewness 
effects become important in the PDFs for the current 
best data sets. 

This paper is organized as follows. In Section [n| we de- 
fine the background spacetime and the precise form and 
distribution of the inhomogeneities, and review the basic 
formalism needed to compute the weak lensing conver- 
gence. In Section [lll| we develop the probabilistic for- 



malism to calculate the lensing magnification PDF. In 



Section IV we derive the analytic approximation for the 
lensing PDF and its mode and dispersion. In Section [V| 
we apply our methods to compute the PDFs for both 
the ACDM model and the EdS model for different types 
of matter distributions, and compare our numerical and 



analytical results. In Section VI we will give our conclu- 
sions. In Appendix |A] we give some analytical results for 
the EdS model and in Appendix [B] we prove an identity 
used in the calculations. 



II. SETUP 

In order to concentrate purely on the effects of weak 
gravitational lensing, we will ignore a possible strong 
backreaction (see Ref. [36^ for a definition and, for ex- 
ample, Ref. [sillMllMlllfflSIlSlliailllia 

[47| SHI [50] for a discussion). In particular we 
will assume that the spacetime of the inhomogeneous 
universe is accurately described by small perturbations 
around the FLRW solution whose energy content and 
spatial curvature are defined as Hubble-volume spatial 
averages over the inhomogeneous universe. Following 
Ref. [36 we will call this Hubble-volume average the 



Global Background Solution (GBS), while the cosmolog- 
ical background solution actually obtained through the 
observations will be called the Phenomenological Back- 
ground Solution (PBS). The inhomogeneities we will in- 
troduce will cause the PBS depart from the GBS, and 
so this work falls into the category of weak backreac- 
tion. Also other phenomena, such as a local Hubble bub- 
ble [2i[271[2ai23l3Ql|3ll|32lE3llMll3^ or redshift 
effects [171 HH [IE [SI] could be studied within the weak 
backreaction scenario. In a companion paper [24 we have 
extended the current work to include a local Hubble bub- 
ble, but here we will ignore this possibility and we also 
neglect all redshift effects [T. In what follows, we will 
first introduce the GBS and the detailed form of the in- 
homogeneities. Then, after defining the parameters that 
describe our model universe, we will briefiy review the 
machinery necessary to calculate weak- lensing effects. 



A. Global background solution 

In agreement with CMB observations we will focus on 
spatially fiat models. Moreover, since we are only in- 
terested on the late evolution of the universe {z < 1.6) 
we neglect radiation retaining only the contributions 
from (dark and baryonic) matter and the dark energy 
in the form of a cosmological constant. The parame- 
ters that specify the GBS will be therefore be 1^m,o and 
Hq. As special cases we will consider in particular the 
ACDM model with 1]m,o = 0.28 and the EdS model with 
^M,o = 1- The evolutions of the Hubble expansion rate 
and of the density parameters as a function of redshift 
are given by 



H{z) 
VLm{z) 



Ho {^Mfi (1 + z f + ^Afi 



1/2 



!^M,0 (1 + zf 

1 - ^m{z) , 



HHz) 



(1) 

(2) 
(3) 



where the subscript will denote the present-day values 
of the quantities throughout this paper and Hq = lOO/i 
km s-^ Mpc-^ Substituting in Eq. ^ H = d{t)/a{t) 
and 1 -\- z = ao/a{t) we obtain the equation we have to 
solve in order to find the time evolution of the GBS. The 
observables we will be interested in are the angular di- 
ameter distance, the luminosity distance and the distance 
modulus: 



7^ / N ^ dz 

""'^'^ = YT-zlm 

Dl{z) = {l + zfDA{z) 



m{z) = 5 log 



10 



Dl{z) 
IMpc 



25. 



(4) 
(5) 
(6) 



Analytical expressions for these quantities in the EdS 
case can be found in Appendix [A] 



3 



B. Inhomogeneities 

We will first describe the statistical distribution of the 
inhomogeneities within the GBS discussed in the previ- 
ous section. The more precise properties of the inho- 
mogeneities, or halos, will be described in the second 
subsection and in the last subsection we summarize the 
parameters that describe our model universe. 



1. General statistical properties 

The inhomogeneity scale introduced by the halos will 
be much smaller than the Hubble radius, dc <C c/Hq^ so 
that the GBS can perturbatively describe the spacetime. 
In particular we will assume that redshifts can be related 
to co-moving distances through the GBS. In our setup the 
halos are randomly distributed with a co-moving number 
density ric and all the matter in the universe is within 



these objects (see Eq. 13). We will later generalize this 



picture for a continuous mass distribution of halos. 

Let us now consider a co-moving volume V much larger 
than the characteristic scale so that the total number 
Nh of halos in V is large, Nh — ricV 1. Let us first 
estimate the average co-moving distance dc between the 
halos. To this end we need the probability P{k; v) of 
having k objects in the co- moving volume v: 



P{k;v) 



Nh- 



k 

k ) \Nh) 

k\ ' 



V \Nh- 

v) 



1 



ricV\ 



Nn-k 



(7) 



where we first wrote the binomial probability for p = v/V 
and then took the limit Nh ^ 1. We have thus found 
that the number of halos in the volume v is distributed as 
a Poisson random variable of parameter UcV which is, as 
expected, the mean number of halos. In particular, mean, 
mode, variance and skewness functions characterizing the 
distribution are: 



k 
k 
a 

71 



lncv\ 



(8) 



where the notation [x\ refers to the floor function cor- 
responding to the largest integer not greater than x and 
7i is the skewness. These are of course well known sta- 
tistical properties of the Poisson distribution. However, 
the skewness of the distribution is so crucial to our re- 
sults that we illustrate this behaviour as a function of ricV 
in Fig. [l] It is evident that the distribution is strongly 
skewed for small UcV^ whereas for increasing UcV the 
mode starts to approach the mean (marked by a vertical 
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FIG. 1: Poisson distribution with maximum normalized to 
unity for four different values of ricV. The vertical dotted 
lines mark the mean. 



dotted line), skewness goes to zero and the distribution 
approaches the normal distribution. In particular the 
difference between mode and mean for low values of ricV 
will be of central relevance for us in what follows. 

Let us now return to the evaluation of dc. Evidently 
UcV is a stable expectation for the number of objects 
in the volume V. To see this imagine that the volume 
V is embedded within an even bigger volume W. The 
probability of having k objects within V is then given 
by the Poisson probability P{k; F), but because ricV 
1 this is is well approximated by a normal distribution 
with mean k = ricV and variance cr = {ucVy^'^. For 
UcV ^ oo we therefore obtain a/k 0; that is, the 
relative fluctuations around the expected value will go 
to zero and we can simply write Nh = UcV. The latter 
property that the Poisson distribution tends to a delta 
function for large parameters will be important when we 
will discuss the lensing bias with respect to the size of 
the data samples of observations. 

Because the objects are randomly distributed, the 
probability w{r)dr that the nearest neighbour to an ob- 
ject is at a radial distance between r and r -\- dr equals 
the probability of having no particle within r times the 
probability of having at least one particle between r and 
r + dr: 

w{r)dr = P{0]v)P{> l]dv) 

= P{0;v){l- P{0;dv)) 

= e-^^^(l - e-^^^^) = Tie dv e-^^^ , (9) 

where v = Airr^/S. Therefore, the average co- moving 
distance dc between nearest neighbours is [57] : 



dc 



r r(4/3) _i/3 



0.55n-^/^ (10) 
where we put Vmax = oo using the fact that V ^ ^c^- 
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2. Detailed description of the halos 

We take our halos to be spherical and characterized by 
a co-moving radius R where the density of the halo goes 
to zero, so that beyond this distance a halo does not 
contribute to the weak lensing convergence. The den- 
sity profile within the radius R can be taken to be any 
smooth function like a gaussian (in which case R 3cr), 
a singular isothermal sphere (SIS) or the Navarro- Frenk- 
White (NFW) profile [58 . Ah matter is taken to be in 
these halos, and so the distribution of matter within the 
co-moving volume V is: 



Pm 



Nh Nh 



(11) 



where pm Pm,o(1 + z)^ with pM,o 3H^Qm,o/{87tG). 
Throughout this paper we use the overbar to denote a 
quantity corresponding to the GBS. The halo profile (p{r) 
must be normalized appropriately to get the average den- 
sity. From Eq. (11) we find: 



Nh 



1 = ^Xev.^"-^/,""^' 

(pdV=^, (12) 



where (pj = (p{\r — rj\) and in the second equality we 
used the fact that Nh ^ 1 in order to neglect boundary 
corrections. Equation (12) ties the halo mass to the co- 



moving halo number density: 



M 



H 



JV 



- 3 



(13) 



and Tic is further related to the average distance dc be- 
tween halos through Eq. (10). In the present analysis Uc 



is a constant with redshift: we will extend this picture 
in Section [III D[ The total energy content of the universe 
and its density contrast can now be expressed as: 



s = 



Pm + Pa 

^ _ Pm - Pm 
P 



(14) 



= ^M 



P 



. Nh 



^M ( X] - l) = ^^E , 



(15) 



where 8m = Pm/pm — 1 and we have defined Sh = 
Qm Sj2i and Se = — 1 = — ^M- The latter gives 
the density contrast in the empty space (voids) between 
the halos, while the former the matter field due to the ha- 
los. The average contrast of a halo is found by averaging 
5h over the volume of a halo: 



(16) 



Eq. (16) will be used to relate the truncation radius R to 



the average contrast at virialization. 

We stress at this point that by "halo" we do not mean 
only gravitationally bound systems, but also non virial- 
ized large-scale structures for which the radius R is not 
rela ted to the virialization contrast of Eq. ([l6| (see Sec- 
tion 



HIE). 



3. Parameters of the model universe 

Summarizing, the parameters that specify our model 
universe are the matter abundance 1^m,o and the Hubble 
expansion rate Hq giving the GBS and the scales R and 
Tie describing the inhomogeneities. Another important 
parameter is the number of supernova observations No 
at a given redshift. As we will see, because of lensing 
effects, the PBS will depend on No and it will reduce 
to the GBS only when the observations cover uniformly 
the entire sky in the limit No oo. The latter limit 
can also be understood as averaging over all the sky. In 



Section III C we will discuss this topic in detail including 
a redshift dependent A^o(^)-function. In Section III A 



we will consider the case where the observations fail to 
cover the entire sky and discuss how these selection ef- 
fects may affect the PBS. Thus, in what follows we will 
study quantitatively the lensing predictions for universes 
described by the parameter sets {flM,Oi Hq, R^ric, No)- 



C. Cumulative gravitational weak lensing 

We will now briefly introduce the tools necessary to 
calculate weak- lensing effects for our setup. For more 
details see Ref. [52j. In the weak- lensing theory the net 
magnification /i produced by a localised density pertur- 
bation is: 



(l-K)- 



(17) 



where k is the lens convergence and 7^ is the shear which 
we assume to be negligible [52 . The shift in the distance 
modulus caused by p then becomes 



Am 



2.51ogio/i - 51ogio(l - 1^) . 



(18) 



The lens convergence k can be computed from the fol- 
lowing integral along the line of sight: 



Jo 



dr 



r{rs - r) V^^ 



(19) 



where is the co-moving position of the source and the 
integral is along an unperturbed light geodesic. The term 
is the Laplacian of the Newtonian potential of the 
perturbation in co-moving coordinates. In the present 
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case then 



V^^ 47rG 



a Sp ■■ 



3 ao 

7:^^M,o — n 

2 a 



Nh 



(20) 

where we have used Eqs. ([l][2| and (15). The sum over 
the halos describes the effect of inhomogeneities, while 
the negative constant of unity can be understood as 
the demagnifying potential of empty space (the "empty 
beam"). Substituting Eq. (20) back to Eq. (19) we get 



nE= j dr G{r, r^) (^^cpj -Ij , 
Jo 



(21) 



where hce the empty beam convergence and k^h is the 
convergence caused by halos, and we have defined the 
auxiliary function G{r^rs): 



G{r,rs) -ftM,o — n 

2 To 



r) ao 



a(t{r)) • 



(22) 



In an exactly homogeneous FLRW model the two con- 
tributions kh and a^^; in Eq. (21) cancel and there is 



no lensing. In our setup, as we will see in detail in the 
next section, the two contributions cancel also if we can 
take the angular average over all sky. In other words, 
given a number of observations Nq^ the GBS is obtained 
in the limit No oo, if no lines of sight are obscured. 
When both these conditions are met, the entire volume 
of the space is seen and the deduced sum ^ ifj averages 
to unity. However, it is evident that a limited size of the 
data set (small Nq) can cause a probabilistic deviation 
of the PBS from the GBS, whereas selection effects can 
lead to systematic deviation of the PBS from the GBS. 
One can imagine several qualitative selection effects that 
could cause an apparent demagnification of the data sam- 
ple, such as rejection of "outliers", intergalactic extinc- 
tion or foreground light contamination. The simplest way 
to estimate such effect is to assume that some fraction of 
the halo mass is "hidden" from the observations. 

Note also that for given parameters (i?, nc), the den- 
sity contrast scales linearly with Qm so that a flat dark- 
energy-dominated model will have smaller lensing effects 
than does a flat matter-dominated EdS model. 



The accuracy of the weak lensing approximation 

Above we defined a cumulative expression for the weak 
lensing convergence hz^ whereas the observable flux is 
actually multiplicatively magnified by the subsequent 
lenses. In practice this works well when the cumulative 
magnification is small. That is, when I75I <C 1 we have 



(23) 



and {k,) = 0, {/a) = 1 and (Am) = are all equiva- 
lent. The accuracy of the weak lensing scheme can be 
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FIG. 2: The re lati ve error Err of the weak-lensing calculation 
defined in Eq. ( 24 ) for the empty beam {a = 0) and for a half 



filled beam {a = 0.5) for both the EdS and the ACDM model. 



tested quantitatively by a comparision with exact solu- 
tions. First, it was shown by Ref. [16 that the weak- 
lensing approximation reproduces the full general rela- 
tivistic results of Ref. |^ to a ~ 5% accuracy along 
a line of sight in an inhomogenous swiss-cheese model^, 
where the difference comes from neglecting higher order 
terms in the expression for hi in Eq. (19). 



Let us next consider the solutions for empty or par- 
tially filled beams with a filling factor a = Pbeam/p- 
Physically this corresponds to the case where a fraction 
1 — a of the mass in the halos is for some reason hidden 
from the observations and the remaining mass fraction a 
is observed as a smooth distribution. Let us define the 
relative error between the weak lensing approximation 
and the exact result for a partially filled beam: 



Err = 



(24) 



In Fig. [2] we show Err in the case a = and a = 0.5 
for the EdS and the ACDM models. The convergence 
factors for the ACDM model were computed integrating 
numerically the results of Ref. [59 . Analytic expressions 
for the exact and weak lensing convergences K,E,a in the 
EdS model are shown in Appendix |Aj Clearly the weak 
lensing and exact results agree to within a few percent 
over the interesting redshift range. In this paper we are 
concerned not with all possible magnifications but with 
the most probable ones. Since these are bound in mag- 
nitude by the empty beam convergence, the results of 
Fig. |2] then suggest that the weak lensing approximation 
should be very good for our purposes. 



1 The authors of Ref. [16] found Am ~ 0.346 at redshift z ~ 1.86 
through a particular grid of LTB-bubbles. They claimed only 
an accuracy of ~ 10% probably because, unlike us, they did not 
have the exact results of Ref. . 



III. PROBABILISTIC STUDY OF WEAK 
LENSING 

Our next task is to compute the probability distribu- 
tion function (PDF) and the most Ukely value of the 
lens convergence hz along arbitrary photon geodesies as a 
function of our model parameters {Qm,Oi Hq, R^ric^ Nq)- 
Since we already know how to calculate remains to 

evaluate the halo-induced part in Eq. (21): 



Jo 



(25) 



where = r{zs). We wish to obtain a probabilistic 
prediction for this quantity along a random line of sight 
to a source located at r^, in a universe with randomly 
distributed halos. One way to solve the problem would 
be to construct explicitly a large enough co-moving vol- 
ume with halos at fixed random locations r^, and then 
compute kh along randomly selected directions in this 
space. Alternatively, one can start from a fixed geodesic 
and construct a random distribution of halos along that 
geodesic. We choose the latter approach, because it al- 
lows us to find fast numerical and analytical solutions for 
the key quantities we are interested in. 

First consider a particular realization 6 of the integral 
in Eq. (25), that is, a particular configuration of the halos 



and a particular line of sight. Because of the finite size 
of the halos, only the A^^i halos with impact parameters 
bj < R from the geodesic contribute to the sum (see 
the sketch in Fig. |3|. Moreover, because the halos are 
small compared with the horizon scale, the function G is 
to a good approximation a constant G{r^rs) ~ G{rj^rs) 
within each halo. Similarly, the time dependence of the 
halo profile will be weak, (p{x^t) ~ (p{x^tj). After a little 
algebra one then finds: 



i^h{0,Zs) ^ ^^Girj.rs) / 



r^^ 2xdx 



^G(r„r,)r(6„t,) 



--(p{x,tj) 



(26) 





FIG. 3: Shown is a co- moving segment of a photon geodesic 
between an observer at O and a source at S. The shaded 
disks represent halos of radius R. 



where "^fjiNi = Nq. Now divide also the impact pa- 
rameter into Nn bins of width A6^, such that we can 
take r(6, t) a constant within each of these bins. In this 
way we can rewrite (27) as 



Ns Nr 

^H{e,z,)^Y.Y.^i 



.G{h.rs)T(hmAri)). (28) 



= 1 m=l 



where hm is the average h within the bin [6^, hjn+i] and 



Ns 



(29) 



i=l m=l 



Now observe that G{r^rs) and T{b^t) are universal func- 
tions for arbitrary values of r and 6, and so all informa- 
tion specific to the particular realization in the equa- 
tion (28) is given by the set of integers {kf^} giving the 
number of halos within the distance and impact parame- 
ter bins characterized by the positions {f^, 6^} and sizes 

Equation ([28| is the starting point of our analysis, be- 
cause it can be easily turned into a probabilistic quan- 
tity. Instead of considering a set of realizations along 
arbitrary lines of sight through a pre-created model uni- 
verse, we can define a statistical distribution of conver- 
gences throug h Eq. (28) because, as we have showed in 



is now characterized by a set of values {rj,6j}, where j 
runs from 1 to A^^i. 

Next divide the geodesic in Ns subintervals of (possi- 
bly variable) length Ar^ R^ such that G(r, r^) can still 
be taken a constant within each interval. In practice one 
can always construct such a division with Ar^ ^ R. 
When these conditions are met, all halos in a given bin 
rj e [r^,r^+i] can be associated with the same distance 



where Rj = R(tj) and tj = t(r,). Our sample realization Section[nBl| the integers i^f^ are distributed as Poisson 



random variables of parameter AA/'^, 



(ATV,^)^- 



(30) 



where ANim is the expected number of halos in the bin 
volume AVim- 



fi to the center of the bin, and (26) becomes 



AN,,, 



UnAVir 



.27TbmAb 



(31) 



Kh{0,Zs) 



Ns Ni 

^G(f„r,)^r(6,,t,), 
i=i 1=1 



(27) 



Equivalently, we can interpret ANim as the mean number 
of collisions of a photon with halos within the z'th r- 
subinterval and within the m'th impact parameter bin. 
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Physically this statistical distribution of convergences is 
equivalent to our original set of realizations 6 averaged 
also over the position of the observer. Thus the statistical 
model explicitly incorporates the Copernican Principle. 

The basic quantity in our probabilistic treatment is 
a configuration of random integers {kim}- The conver- 
gence corresponding to a given configuration is given by 
an equation analogous to Eq. ( 28 ) : 



where Lv,'s are random variables drawn from the Pois- 



Ns Nr 

EE---- 

i=l m=l 



where hziim{zs) = G{ri,rs)T{bm^U) is the convergence 
due to a single halo at distance and impact parameter 
bm- The probability that such configuration occurrs is 
just 



Ns Nr 



= n n = n^'=- 



(33) 



= 1 m=0 



It is easy to see that this probability is normalized to one: 



E ^{fc™} 



n 



im kim = ^ 



(34) 



The expectation value for the convergence is given by the 
probability weighted sum over all possible {kim} config- 
urations. Using the above results it becomes: 



{kh{Zs)) = ^ P{kim}'^H{{kim},Zs) 

^ ^ {kim) f^limj^s) 
im 

= ^ ANiml^iim{Zs) 
im 

Ns 

= G{ri,rs)Ari = -ke ' 



(35) 



That is, the total expected convergence (k) = (i^h) -\- i^e 
vanishes bringing back the GBS result, consistent with 
photon conservation in weak lensing. In going from the 
third to the last line we used the identity 



Nr 

E 

m=l 



27TbmAbm ncT{bm,ti) 
R 



I 



27rnc / dbbr{b,t) = 1, (36) 



which, when calculated in the continuum limit, holds true 
for any functional form of the halo profile (f with the nor- 
malization of Eq. (12) as shown in Appendix p| Finally 



son distribution (30): kim = kim[ANim]' Eq. (37) makes 



explicit that the expected convergence vanishes: the ex- 
pected value of each term of the summation is indeed 
zero. 

Note that in the analysis above, we did not give any 
quantitative criterion for binning the variables r and b. 
The fact that both Ar^ and Abi formally vanish in the fi- 
nal expression in Eq. ( 35 ) already suggests that the exact 



way the binning is done is not important. This is indeed 
so, and one can even show formally that the binning is not 
important as long as it is sufficiently fine to give a good 
approximation for the functions G(r, r^) and r(6, t). One 
can also easily test the effect of binning directly by com- 
paring the final results obtained with different size grids. 
In practice we found that a grid with 10 — 20 points in 
each variable already gives very accurate results. 

The final convergence PDF can be formally written as 

P^i{k,Zs)= ^ P{/e,^}J(/^ - /^({/Ci^},^^)) . (38) 

{kim} 

In the continuum limit this becomes formally a functional 
integral over random integer- valued functions k{r^b). It 
is clear that the most likely configuration which max- 
imises the probability function P{ki^} corresponds to the 
mode: kim L^^^mJ- Moreover, for large ANi^ this 
most likely configuration approaches the mean kim 
ANim- When this is the case the expectation value for 
total convergence vanishes even for a single observation. 
This signals the fact that the PDF then approaches a 
gaussian with a vanishing skewness. 

One can also create the probability distribution di- 
rectly for the shift in distance modulus Am: 



P^\{Am,Zs) 



E ^(^. 



J (5(Am - /\m{{kim}, Zg)) , 

(39) 

where one uses equation (18) to compute Am{{kim}^ Zg) 
from the converegence. Similarly, one can also define the 
PDF form the magnification /a. 



It is straightforward to compute the distributions ([38|) 
and 



39[ ) through a numerical simulation, simply by cre- 
ating a sufficiently large set of random configurations 
{kim} drawn from the probability distribution P|/g.^}. 



Note that the distributions (38) and (39) are formally 



discrete but very dense sets of distributions whose inte- 
gral over hz or Am is normalized to unity. So the quantity 
that one actually computes through a simulation is: 



nAml+Abir^/2 

P^\{Ami,Zs) = / dyP^\{y, 

J Ami-Abiri/2 



zs) , (40) 



note that, using the identity (|36|), we can write the total ^here A^n are some suitably chosen bin widths. The 



convergence for an observation corresponding to a con- 
figuration {kim} simply as: 



PDF (40) is one of the main results of this paper. We 



Ns Nr 



f^{{kim}^ Zg) — ^ ^ ^ ^ f^lim{Zs)(^im ^^im^^ (37) 



i—1 m=l 



wish to stress the simplicity of this result. In order to 
specify the model completely one needs only Eqs. (22) 
and (26) for the functions G(r, r^) and T{b^t) respec- 
tively. After this, P(Am/,Zs) is found from the random 
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sample of magnifications Am computed through equa- 
tions (Isrl) and (IT8|). 



the simple selection effect discussed in Section [II C| To 



see this, rewrite Eq. (45) as 



A. Selection effects 



Ns 



{k{zs)) = ^{a, - l)G(f„ r,)Ar, , (46) 



As was discussed in Section II C , several selection ef- 



fects might bias the convergence distribution, favouring 
overall demagnification. It is straightforward to general- 
ize our probabilistic approach to include many such ef- 
fects. For example, all effects that would lead to a rejec- 
tion of a potential SN observation can easily be described 
by an additional survival probability function. That is, we 
replace 



peff 



(41) 



in our master formula for the magnification PDF (40). 



Here K is a normalization constant that makes sure that 
the final PDF is normalized to unity. The most generic 
form of the survival probability function is 



\Kzm f X J. ^ tin / 



(42) 



This allows the surival function depend on the arbitrary 
local properties along the photon geodesic. Given the 
form (|42|) the proper normalization is easily seen to be: 



3eff 

'{kim } 



=n 



k- ' 



kirr 



where 



(43) 



(44) 



The correct expression for the convergence is still given 
by Eq. (37), with the important difference that the ran- 
dom integers kim are now drawn from the Poisson dis- 
tribution Eq. (|43| with the effective expected number of 

halos: A^^m = ^im 

[ATVfff]. Now {hm) = ATVfff so that 
the average convergence over a large number of observa- 
tions becomes 



(45) 



As expected, the selection biases survive in the data even 
after averaging over many observations. Note that in the 
above analysis we implicitly assumed that the survival 
probability for a light ray going through empty space 
equals unity. 

The actual form of the survival function depends 
on detailed input both from the astrophysical properties 
of the intervening matter distributions (estimate for ex- 
tinction, foreground light contamination, etc.) and from 
the observational apparatus {e.g. detection efficiency of 
SNe-triggering telescopes). A quantitative analysis of 
these issues is beyond the scope of the present paper, 
and we merely show how the method can be applied to 



where 



ai 



Tie 27rbmAbm T{bm, ti) PI 



sur 
im 



(47) 



is the fraction of the mass in the halos that is ac- 
cessible for observations at n. Note that a constant 



psur 

im 



thanks to the identity (36), is a special so- 



lution to Eq. ( |47[ ). If P^^ did not depend on r^, then 
this further reduces to 



(hi) = {1- a) ke 



(48) 



with ai = a for all i. As outlined in Section |II C| and 
showed in detail in Appendix [Aj a can be interpreted as 
the filling factor a = Pbeam/ p of a partial filled beam. We 
see therefore that the results of Ref. [4 can be obtained 
as a limiting case of our approach. 

Other types of selection biases, such as due to an error 
in the estimate for reddening, would not affect the prob- 
ability distributions, but the evaluation of the effective 
magnification factor itself. Also these corrections could 
depend on the mass distribution along the photon path, 
and so, to be as general as possible, one should replace 



Am{{kim}, Zs) + SsAm{{kim}, Zs) 
Amef^{{kim},Zs) (49) 



Am 



in the integrand of Eq. ( |4Q| ). Again, a quantitative esti- 
mate of the size of the bias factor ^gAm is beyond the 
scope of this paper. 



B. Sources with intrinsic luminosity dispersion 

Our results can be easily generalized for an arbitrary 
initial flux or magnitude spectrum of an imperfect stan- 
dard candle. If we know that the source magnitudes can 
be described by a function Pin(Amo), then the observed 
magnification PDF is obtained by the convolution of the 
initial PDF and the weak lensing PDF 

P{Am,Zs) = J d^Pwi(y)Pin(Am-^) 

= Yl P{k,m}Pin{Am-Am{{kim},Zs)). (50) 

{k.m} 

Alternative expressions where one or both distributions 
are replaced by the probability distribution in the flux 
are trivially obtained by a simple change of variables. 
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C. PDF for a sample of No observations 



D. Arbitrary mass-distribution of halos 



The effective distribution for a binned sample of No 
supernova observations can formally be expressed as an 
iterative convolution [12] : 

PNoizs^Am) = J dyPNo-ii^s^y) X 

xP^{NoAm-{No-l)y). (51) 

Here Pi is the normalized fundamental PDF that may 
include selection effects and a convolution over the ini- 
tial distribution. Alternatively, one can compute Pnq 
directly from the fundamental PDF by creating a large 
number of sets of No random realizations from it, and 
creating a normalized distribution for the average Am 
within these sets. 

A third way is to create the Pnq distribution directly 
in the initial simulation, bypassing the calculation of the 
fundamental Pi altogether. This reduces computational 
time and improves accuracy, of course with the price of 
limiting the amount of information available. If we label 
the configurations within a given sample of observations 
by 5, the mean convergence after No observations is: 



f^No ({^im}) 



No 



Ns Nr 



i=l TO=1 ^ 

Ns Nr 

EE^i-(%^-^^*-)' (52) 



= 1 m=l 



where, given that each of the No observations are inde- 
pendent, we have used the fact that independent Poisson 
variables with the same weight sum exactly into a Poisson 
variable of parameter given by the sum of the individual 
parameters: kim,No = kim[NoANim]- Pno is then given 
by Eq. ([39| where the magnification for a given configura- 
tion {kj^ is computed from the convergence of Eq. (52). 



We followed this last approach when creating the illus- 
trations relevant for the Union Catalog and JDEM data 
in Section [Yl However, note that including the selection 
effects in general does not commute with taking the av- 
erage over the observations. In other words, if the No 



measurements are correlated then we cannot use Eq. (52 ), 



but we have to start from the fundamental PDF. On the 



other hand, Eq. (52) displays explicitly the effect of the 
size of the data sample: even when has a skewed PDF 
and a nonzero mode for No = 1, for large No the distri- 
bution approaches a gaussian and eventually converges 
to a ^-function at zero convergence, as it is clear from 
the properties of Poisson variables discussed in Section 

HlBTl 



We have so far expressed all our derivations assuming 
there is only one halo type. At this point it is straight- 
forward to generalize all our formulas to the case where 
one has a continuous distribution of halos. Let us as- 
sume that the halo mass distribution is given by some 
dimensionless function /(M, z), which is related to the 
co-moving number density n(M, z) through the standard 
relation 



n(M, z) 



Pm,o clq 
M 



f{M,z). 



(53) 



We also assume that /(M, z) is normalized to one, say: 



dMi5/(Mi5,^) = 1, 



(54) 



where we have introduced the usual dimensionless mass 
parameter M15 = M/{10^^h-^MQ). That is, our def- 
inition for /(M, z) is analogous to the usual Press- 
Schecter [60 and the Sheth-Tormen [61 mass functions. 
However, our /(M, z) is not expected to be quantitatively 
similar to the PS- and ST-functions, since the latter are 
designed to model the density of virialized systems, while 
we are also interested in effectively describing other large- 
scale structures (see the next Section HIE). 

The first step to extend our treatment is to generalize 
the normalization of the halo profile cp: 



4:71 



I 



R{t,M) 



ip{b,t,M) b^db 



f{M,t) 



(55) 



n(M,t) ' 

where instead of the redshift z we have used the time 
along the geodesic, /(M, t) = f{M,z{t)). The halo mass 



follows as in Eq. (26): 



is now given by Eq. (53) and the definition of r(6, M, t) 



r{b,t,M) 



I 



R{t,M) 



2xdx 



62 



ip{x,t,M). (56) 



Finally we discretize the function /(M, z) in the same 
way we discretized the impact parameter and the co- 
moving distance and the final result is: 



Ns Nr Nm 

^({^imn}? ^s) ~ ^ ^ ^ ^ ^ ^ ^limn(^s) 
i=l m=l n—1 



X ( kir. 



(57) 



where Kuran{zs) = G{fi,rs)T(bra,ti,Mn) and Mn is the 
mean M in the n'th mass bin. The parameter of the 
Poisson variable kimn = kimn[ANimn] now is 



(58) 



The probability distribution for the shift in distance mod- 
ulus Am now becomes 

x6{Am- Am{{kimn},Zs)) . (59) 



10 



One can similarly extend the systematic effects to be de- 
pendent on the mass, such that the survival probability 



whereby AiVf^ ^ ^Nf^n in Eq. (|44|), 



sirable to improve the clustering algorighm, and we can 
imagine several ways to do it within our approach. First, 
our most general result in the form of Eq. (57) allows 



and also Am Ameff({/^imn}7 ^s) in Eq. (49). Finally, 
Eq. (40) formally retains its present form. Eqs. ( 57p9 ) 
can be used to compute all the interesting quantities for 
an arbitrary halo mass function /(M, z). However, for 
simplicty, we will from now on concentrate only on single- 
mass functions in this paper. See the model of Ref. [24] 
for an example with two types of halos. 



E. Discussion 

Our stochastic analysis of lensing presented above is 
based on using an unclustered Poisson distribution of 
halos. Let us now discuss the applicability of this as- 
sumption and the possible ways to improve the clustering 
algorithm. 

Ref. [6] states that the lensing properties of a universe 
made of point masses are independent of their masses and 
clustering. The statement was confirmed numerically in 
two steps. First, it was shown that random distribu- 
tions of point masses of M©, lO^^M© and lO^^M© share 
the same lensing properties. Second, Ref. [6 compared 
a universe made of point mass galaxies with a universe 
made of uniform density galaxies composed of point mass 
stars. The authors concluded with their Fig. 10 that, for 
an intergalactic separation of 2 Mpc (which corresponds 
to a mass of ~ 10^^ M©) and a radius for the uniform 
density galaxies of 200 kpc, the above two universes have 
the same lensing properties. 

These results, however, do not apply for extended ha- 
los, and anyway do not prove that clustering can be ne- 
glected for objects with masses larger than ~ lO^^M©. 
On the contrary, we will show in Section |V| that the lens- 
ing PDF changes significantly when one considers halos of 
the scale of large superclusters of galaxies (compare Fig.js] 
with Fi g. [9|). Let us stress again, as we pointed out in 
Section IIB 2[ that by "halo" we do not mean only gravi- 
tationally bound systems, but (actually the surface mass 
projection of) any unvirialized sufficiently localized large- 
scale structures. The idea is that an unclustered Poisson 
distribution of large-scale structures can give a qualita- 
tive estimate of the lensing effects produced by the actual 
clustering of galaxies into filaments and walls. This idea 
was carried out further in the companion work of Ref. [24 
where structures of average separation of dc = 100 
Mpc (which corresponds to a mass of 6- lO^^/i^^M©) and 
radius of Rp = 10 Mpc were considered. 

Of course randomly placed spherical structures is but 
a crude approximation for the actual 3-dimensional web- 
like structure of clusters and filaments. However, all that 
matters for weak lensing are the projected 2D-surface 
mass densities over a series of spatial slices from the ob- 
server to the source. For these projections the difference 
between a realistic 3D-structure and the crude meatball 
model is less distinctive. Nevertheless, it would be de- 



us to deal with a generic halo mass distribution function 
/(M, z). In the most conservative approach one would 
identify / with some known clustering function, such as 
the Press-Sechter distribution. However, to incorporate 
also the noncoherent structures one could use instead an 
effective form, say /(M, ~ N^fps{£,M,z)^ where 
is a normalization factor. For ^ < 1 such form would 
in a simple way model the merging of smaller halos into 
larger systems of (unvirialized) clusters. Another, and 
probably a more accurate method to describe very large 
scale clustering would be to use a master probability dis- 
tribution function to modulate the average matter back- 
ground density on the scale O(IOO) Mpc, along the lines 
of sight, together with a normal local clustering fuction 
/(M, (5p(z)), where dp is the modulation around the 
GBS density. This method also allows an easy way to 
implement a local void around the observer. Yet an- 
other way to improve the clustering algorithm would be 
to introduce non-spherical structures, such as cylinders to 
better describe filaments and walls. However, we suspect 
that while such extension would apparently help gaining 
a more realistic 3D-distribution, it would not be crucial 
for the modelling of the surface density projections. 



IV. ANALYTIC RESULTS FOR Pwi(Am, Zs) 

In this chapter we will derive an analytic approxima- 
tion for the magnification PDF defined in the previous 
Section, in particular we will focus on mode and vari- 
ance. A reader not interested in details can skip to the 
final results presented in Eqs. (83), (87) and (94); the 



validity of these expressions will be tested numerically in 
the Section |V| The procedure will consist of two steps of 
resummation over the variables bm and r^. Let us first 
introduce the central idea for these summations using a 
generic example depending on one set of variables. 



A. Linear combination of Poisson random variables 

Consider a collection of independent Poisson random 
variables ki with parameters A^, ki[Xi]. As we have seen 
in the derivation of Eq. (52), the sum k = ki is then 
a Poisson random variable with parameter A = J^i^i- 
However, we wish to study the distribution of the quan- 
tity 



S — y ^ Zi ^^^[A^] , 



(60) 



where the weights Zi are positive real numbers. When 
the weights Zi are different, no simple exact expression 
exists for the probability distribution of S. We shall 
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adapt the following approximative procedure: 

i i 
i i 

= Zk[X] =S, (61) 

where we have defined Z = max{Z^} and Zi = ZijZ. The 
new variable S is Poisson distributed with the parameter 
A = X1^2;^A^. The idea of the approximation made in 
the third step of Eq. (61 ) is to use the normalized weight 
distribution zi to define a new set of random variables 
Zi\i that favour the terms that give the largest contribu- 
tion to the original sum S. In this way, the PDF for S 
should provide a reasonable approximation for the mode 
and the skewness of the actual distribution. It easy to 
show that this procedure preserves the mean and works 
exactly for the trivial case of constant Zi. However, the 
approximation distorts the variance of the distribution: 



(62) 
(63) 



Thus the dispersion a for S should be be corrected with: 



B. Summation in hr, 



(64) 



Let us now apply the results of the previous Section to 
reduce Eq. (37). Here we have two independent indices to 



be accounted for, and we shall resum them sequentially 
starting from the summation in hm- Writing only the 
relevant pieces we have 



where = max{r(6^, t^)} and ki is the new effective 
Poisson distributed random variable with an effective pa- 
rameter AA^^: 

Ar Ar 

^^i = Y^Yl ^^^ruAhm ric T{hm, U) = ^ , (66) 



where we again used Eq. (36). Inserting the approxima- 
tion (65) in Eq. (37) we then find 



Ns 



A. = ^G(r„r,)(f,fc,[AiV,]-Ar,) . 



(67) 



Before resumming the r-variable we pause to give esti- 
mates for the quantity f . These obviously depend on the 



choice of the halo density profile. First consider the uni- 
form density halo with (/:?uni = 3/(47r i?^nc). This implies 



Jb 



^ 2xdx 
b \J — IP' 



3VR^ - 

27TR^nr 



so that 



Tuni — runi(07-^) — 



(68) 



(69) 



Second, for a correctly normalized gaussian 



one finds 



- gau 



2nc 71 



These results suggest we parametrize our as: 
T-^=n,^R^^Ql, 



(71) 



(72) 



where Ri = R{ti). That is, F^ correponds to an effective 
mean free path of a photon at r ~ r^. The effective 
distribution parameter AA^^ now becomes: 



ANi = nc7TR^,Ql An 



(73) 



For example for the gaussian profile above it is Q^p = 
V2/3 0^ OAT. Equations ^ and ^ specify the dis- 



tribution (67) completely in terms of a single effective 



parameter Q^^ which depends on the particular density 
profile (f chosen. This method of estimating analytically 
fails for profiles that are singular at the origin, such 
as the sigular isothermal sphere (SIS). In these cases Q^^ 
can be fitted from the numerical distribution. In any 
case, the essential point is that Q^p can be taken as a 
constant, and the general rule is that the more peaked 
the profile is, the smaller effective one finds. Table [l| 
at the end of this Section shows the numerical values used 
in Section [Vl 



C. Summation in n 

We now repeat the resummmation approximation for 
Eq. (67). Again writing only the relevant term, we get 



Ns 



= ^^^(r.) k,[AN,] ^ H{rs)k[Nc] , (74) 



where Hi{rs) = G(ri,rs)F^ and again H{rs) = maxji^^}. 
Finally 



Ns 

Nc = Yh,AN,, 

i=l 



(75) 
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with hi = Hi/H. Inserting the expression Eq. (66) for 
ANi into Eq. ([75]) we get 



Ns 



(76) 



where we again defined G{rs) = max{G(ri, rg)} and 
gi{rs) = Giri^Vs) /Givs)' The function g depends very 
weakly on Vg and one finds that numericahy 



E 



dx g{x) 



Finahy one can show that to a very good accuracy 
G{rs) 



(77) 



(78) 



where R = R{t{rs/2)). Inserting these results back to 
Eq. ( 75 ) we finally get the estimation: 



Nc^ric-TsE- ttR^QI . (79) 

Physically Nc corresponds to the expected number of 
collisions with halos in a tube of effective radius RQ^^ 
and effective length rsE connecting observer to source. 
Our final result is an effective convergence function which 
is Poisson distributed with parameter Nc' 



n{k,Zs) = 1 



k[Nc 
Nc 



(80) 



Note that the average convergence for this distribution 
vanishes as it must. 



D. Final analytic result 



Both our exact result of Eq. ( 37 ) and our approximate 



result of Eq. ( 80 ) are representative of the actual PDF 



for magnification and convergence. Often one is rather 
interested in the probabilistic interpretation of a set of 
No identical, or sufficiently similar observations. For ex- 
ample, one typically introduces some binning of the data 
points, and one would like to know what is the most likely 
value and the dispersion for such an effective observable. 
Repeating the derivation of the result of Eq. (52), we can 



extend Eq. ( 80 ) to a sample of Nq observation by replac- 



ing Nc with the total effective number of collisions: 

Nt{z) = No{z)Nc{z) , (81) 
SO that the convergence for a configuration k is 



<k,Zs)=^E (l-^) • 



(82) 



Finally, the approximate probability distribution func- 
tion in magnitudes is: 



P„,(Am,z,) = e-^--^, 



(83) 



where k solves Am = b\og^Q{l — n{k^ Zg)). The dispersion 
of Eq. ([83| has to be corrected by the factor oj to be 
derived below in Section [IV F[ 



The distribution of Eq. (83), similarly to the exact re- 



sult of Eq. (39), is discrete. For the latter the chosen bin 
width of Eq. (40) is irrelevant because the PDF is a 



very dense in its domain. This is not, however, the case 
for Eq. ( 83 ) and so we have to give an explicit prescription 



for the binning. In order to estimate the order of magni- 
tude of Ai)in^ we evaluate the convergence caused by one 
halo placed at half way between observer and source and 
hit with impact parameter RQf. 



Ri^G{rs/2,rs)T{RQf,t{rs/2)). 



(84) 



The bin width will then be A 



bin 



51ogio(l - i^i) 



E. Mode 

The mode for the effective convergence is obviously 

[Nt{z)\ 



kb{z) = ke{z) ( 1 



Nt{z) 



(85) 



We denoted the mode by hzsiz) to emphasize its mean- 
ing as the bias away from the homogenous limit = 0. 
The final step in our derivation is to remove the discrete- 
ness of the probability distribution by introducing the 
continuum approximation x(^) fo^" the fioor function: 



2x 



< X < 1/2 
X > 1/2 



(86) 



which is obtained averaging upper and lower boundaries 
of [xj/x. The final expression for the bias is then given 
by: 



KBiz) = KE{z) (l-xiNriz))). 



(87) 



We shall see in Section [V| below that this simple analytic 
approximation can reproduce the mode of the numeri- 
cally simulated magnification PDF ( [4Q| ) for different red- 
shifts, for different halo profiles and for different values 
of No- Finally, if we include the probabilistic selection 
effects through the parameter a (see Sec. Ill A), we get 
simply 



When Nq is very large, x{No) 
as expected from Eq. (48). 



axiNriz)) 
1 and Kb 



F. Dispersion 

In order to estimate the dispersion of the magnifica- 
tion PDF we first need to evaluate the streching factor 
Eq. (64) for our setup. We define 



(89) 
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where the factor cOr comes from the r^-summation and the 
factor uJb from the 6m-summation. uor is roughly given by 



dxg{x)h{x) 
dxg{x) 



0.77. 



(90) 



Like the factor E above, cOr is an almost universal con- 
stant, essentially set by the form of the function G(r, r^). 

depends on the profile (p and is more difficult to es- 
timate. We find that it is reasonably well approximated 
by the ratio: 



T{RQ^,t{rs/2)) 
mrs/2) 



(91) 



For example for the gaussian profile one finds a;^ ~ 
so that CO 0.53. Similarly to the evaluation of Q^^ 
our derivation does not apply for singular profiles and 
the value of oo has to be obtained from the numerical 
distribution. The essential point is again that co can be 
taken a constant. 

Now consider the limit Nt ^ 1, so that the Poisson 



distribution ( 83 ) can be approximated with a normal dis- 
tribution. Taking into account the stretching factor co we 
have (up to normalization): 



Pwi(/c) cx exp 



2^ 



(92) 



where = ^/Nrix). The corresponding dispersion in 
magnitudes is found by evaluating Eq. ( 82 ) at /c = Nt ± 



2.5 log 



.10 



^/NT{z) - KE{z)iJ 
5 I^e{z)uJ C^Am,leff 



(93) 



We see that the dispersion scales as expected with Nq- 
This derivation makes sense for aAm,No ^^^^ when 
^c{z) <C 1 given enough observations such that Nt ^ 1 
holds. This is actually the way Ref. [12] defined the ef- 
fective dispersion for a singular observation, so that 



5 k.e{z)u 



(94) 



From Eq. (79) it is clear that crAm,ieff depends on the 
ratio R^^ = co/Q^ and not on and co individually. 
Summarizing, is about the bias and about the 
dispersion, two independent degrees of freedom. In Ta- 
ble [J we list the values we will use in the next Section. 



We wish to stress that Eq. (94) and Eq. (87) depend 
analytically upon the parameters of the model universe 
and provide an easy way to estimate magnification bias 
and dispersion. 



TABLE I: 


Effective parameters Qc^ 


and R^. 


Halo profile cp 




Rip 


Uniform 


0.90 


0.9 


Gaussian 


0.53 


1.3 
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FIG. 4: On the top, binned number of observed SNe for = 
0.1 for the Union Compilation of Ref. 3 of 307 SNe. On the 
bottom, binned number of observed SNe for = 0.1 for a 
JDEM-like survey of 2000 SNe, see Ref. [2]. 



Finally, note that the PDF of a single supernova is 
not a gaussian of standard deviation crAm,ieff and there- 
fore we will introduce in the next Section the full width 
at half maximum (FWHM) as an indicator of the dis- 
persion for skewed distributions. The applicability of 
Eq. (94) is indeed restricted only to data sets for which 
the condition Nt{z) = No{z)Nc{z) > 1 holds. In the 
upper panel of Fig. [4] we show the binned data for Nq 
from the currently largest SNe-sample, the Union Com- 
pilation of Ref. [3 . The data set is still quite sparse 
for z > 1, which suggests using the approximation (94) 



might not be applicable there to model large structures 
with Nc{z) < 1. The lower panel in Fig. H] shows the 
simulated data for the future JDEM survey. As we shall 
see in Section |V| almost no bias remains in the effective 
PDF for the JDEM-dataset, even within a universe with 
very large structures. 
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V. RESULTS 

Here we will present numerical and analytical results 
for the magnification PDF. We will first consider a model 
similar to the one studied by Ref. [12] in order to compare 
our method with their ray-tracing simulations. In the 
next subsecion we will consider the possible effects of 
large-scale clustering on the magnification PDF. 

Our numerical results are computed with turboGL, a 
very fast Mathematica code based on the formalism de- 
veloped in Sect. |III[ and publicly available: on an ordi- 
nary desktop/laptop computer it takes less than a second 
to compute a PDF with the usual statistics of 10^ con- 
figurations. This has to be compared to expensive ray- 
tracing techniques. Thanks to this performance we can 
use much better statistics in our analysis and simulations: 
10^ configurations already give rise to very smooth his- 
tograms and curves. We would like to stress that since 
turboGL is so fast, the analytic results of the previous 
section are not used as a tool for the analysis, but rather 
to provide analytic insights into the problem. 
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FIG. 5: Shown are the magnification PDFs for a source at 
z — 1.^ for the ACDM model of Table|ll]with SIS and gaussian 
halo profiles. The SIS case has been evaluated for both No — 
1 and No — 50 SNe measurements. Each histogram has a 
statistics of 10^ realizations. Also shown (filled circles) are the 
approximate PDFs given by the distribution of Eq. (92). For 
the PDFs relative to No — 1 only the mode of the distribution 
is shown. 



A. Comparison with Holz &; Linder (2005) 

Holz and Linder [12 computed lensing effects in a 
ACDM universe where all matter was assumed to exist in 
the form of a random distribution of spherical halos the 
size of large galaxies. Ref. [12 used the ray-tracing simu- 
lation developed in Ref. [6 from which we also obtained 
the value of the parameters Rp and dc and the density 
profile used^. Rp is the proper radius which is assumed to 
be constant, while the value of dc fixes, through Eq. (13), 



the halo mass. The value of the parameters used in this 
Section are summarized in Table HIl 



TABLE II: Setup to be compared to the results of Ref. [12], 



Quantity 



Value 



h 

dc 
Rp 
Mh 

Halo profile 



0.28 
0.7 

2 Mpc 
200 kpc 
1.8 • lO^^Mf^ 
SIS 







In Fig. [5] we plot the magnification PDF for a source 
{No = 1) at redshift z = 1.5 in the case of SIS and gaus- 
sian halo profiles. Both PDFs are visibly skewed, but the 
bias is still realatively weak, thanks to the high degree 
of homogeneity given by the parameters of Table [llj The 
SIS-PDF has both a larger bias and a more pronounced 



tail at high-magnifications than does the gaussian PDF. 
This was expected because the strong cusp of the SIS 
leads to more high-magnification events compensated by 
a drop at small magnifications. Also shown in Fig.[5]is the 
SIS-PDF relative to a set of No = 50 SNe measurements: 
as expected skewness and dispersion are reduced show- 
ing that, as discussed in Section |III C[ the distribution 
approaches a gaussian for large A/'o, eventually reducing 
to a J-function at zero convergence. These results are in 
good agreement with Fig. 1 of Ref. [12], which shows that 
our statistical weak lensing approximation can well repro- 
duce the results of ray-tracing calculations. In Fig. [5] we 
also show (filled circles) the analytic effective probability 



distribution developed in Section [LVj As was discussed in 
subsection |IV F[ the effective PDF distribution is a use- 
ful quantity only for Nt ^ 1, and so we plotted the full 
effective PDF only for the SIS case with No = 50, where 
Nt — 37. The agreement with the exact results is re- 
markably good. For the No = 1 cases the point density 
is not sufficient to model the full PDF, but the modes of 
the distributions are still well produced by the analytic 
approximation (single filled circles on No = 1 curves). 

Besides the standard deviation and similarly to 
Ref. [12 , we will use the full width at half maximum 
(FWHM) as an alternative indicator for the dispersion 
of a skewed distribution. In particular we define three 
different indicators: 

FWHM_ 



1.18 

FWHM+ 



The values of the parameters used by Ref. 12 were not stated. 
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2.35 



(95) 



15 



0.06 

0.05 

0.04 

^ 0.03 
E 

0.02 
0.01 




position of 
tiie mode 




FIG. 6: Shown are the bias (top) and the dispersion (bottom) 
of the magnification PDF in the ACDM model of Table [HI as 
a function of the redshift of the source for No — 1- Filled 
circles mark the values obtained from the numerical simula- 
tion (omitted for clarity in the V- and V-^ curves). crAm,ieff 
is calculated from Eq. (94). 



where FWHM_/^ are defined as the distances from the 
left/right edges to the mode of the total FWHM and so 
FWHM_ + FWHM+ = FWHM (see Fig. [lO] for an illus- 
tration). We found these indicators particularly useful 
because for a gaussian distribution one has V- = V-^ = 

V = a where a is the standard deviation and so depar- 
tures from this limit signal the presence of skewness in the 
distribution. Mode together with will describe the 
skewed peak of a PDF, while v and a its dispersion. The 
idea is that for skewed distributions and small datasets 

V is more meaningful than a, which is affected by the 
long high-magnification tail which is of little importance 
if we deal with small data sets. This is actually the rea- 
son why we focused on the FWHM to characterize the 
skewness: the third standardized moment (the skewness 
7i), for example, would have been again sensitive to the 
long and low high-magnification tail. 

In Fig. [6] we show the bias and the dispersion for the 
SIS profile for Nq = 1. This figure can be compared with 
Figs. 5 and 7 of Ref. [12 . While the bias introduced by 
the skewness of the distribution is moderate as shown 
from the top panel of Fig. |6j there is a clear deviation 
from gaussianity, as one finds that a ^ 2v. Moreover 
c^Am,ieff very well reproduces the numerical results for a 
and consequently overestimates the dispersion v. 

Fig. [7| is similar to Fig. |6] but now constructed for 
No = 50. The bias introduced by the skewness of the 
distribution is now greatly reduced. Moreover the dis- 
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FIG. 7: Same as Fig. [6] but for No = 50. The dispersions v 
and a are now closer than in Fig.[6]showing that, by increasing 
the number of observations, the lensing PDF approaches a 
gaussian. 



persions v and a are converging (<j ~ 1.3 v) showing that 
averaging over a large sample of measurements can re- 
duce the effects of the skewness of the distribution. 

In order to investigate the convergence of the lensing 
PDF to a gaussian, we studied in Fig. |8] the dependence 
of bias and dispersion upon the data set size No- The 
bias (top panel of Figj8| is plotted together with the 
analytic result of Eq. ^STf: k,b scales like 1/No and so we 
conclude that the actual bias converges slower than that, 
in agreement with Fig. 2 of Ref. [12 . The bottom panel of 
Fig. [8] instead shows that a is well described by crAm,ieff, 
which means that a scales, as expected, like 1 / V No • 
This plot agrees again with Fig. 4 of Ref. [12^ and shows 
that a converges to v only for very large No- We stress 
at this point that this strong non-guassian signature is 
not a general feature but is due to the strong cusp of 
the SIS profile and is largely reduced with other density 
profiles like the NFW or the gaussian. 

As a technical note, we point out that to obtain the 
non-gaussian signature of Fig. |8) we had to increase the 
binning in the impact parameter in order to accurately 
resolve the cusp. This of course increases computational 
time and so it is worth asking what is the impact of a less 
accurate binning. Our conclusions are that losing that 
non-gaussian signature is not an important problem be- 
cause the skewness of the distribution is underestimdited 
only when the latter is already very close to a gaussian 
and thus at worst leads to a slight conservative underes- 
timation of the bias for large data sets and SIS profiles. 

The results of this Section show that our statistical 



16 



CO 

^ 0.01 



0.1 - 



CO 

E 0.04 



0.02 



z = 1.5 



position of tine mode 




TABLE III: Parameters modelling large-scale structures. 



FIG. 8: Shown are the bias (top) and the dispersion (bottom) 
of the magnification PDF in the ACDM model of Table [ll] as a 
function of the size No of the data set for a source at redshift 
z — 1.5. Filled circles mark the values obtained from the 
numerical simulation. crAm,ieff is calculated from Eq. (94) 



and Kb from Eq. (87) 



model, based on cumulative weak lensing, very well de- 
scribes the most important effects of lensing by the in- 
homogenous matter distributions. However, the model 
discussed in this Section, defined by the parameters in 
Table [llj neglects all effects of large scale clustering on 
the weak lensing. Such clustering is nevertheless clearly 
present in the form of large voids, superclusters and fil- 
amentary structures seen both in the large scale simula- 
tions and in the actual redsift survey datasets. We will 
next try to get an idea of the potential impact of these 
effects by considering a universe made out of much larger 
halos. 



B. Large-scale structures 

In this Section we study the lensing effects caused by 
very large-scale structures. For simplicity we still consid- 
ered universe models made of a statistical distribution of 
a single type of halos. We performed several simulations 
with different halo masses and found that the bias effects 
start to become significant when the halo masses are at 
least of the order of the largest gravitationally bound 
superclusters: M > lO^^/i^^M©. However, even larger 
nonvirialized structures can exist and to illustrate their 
effects we chose models with even larger masses. The de- 



tails of the models are given in Table HI We will consider 
both the ACDM model with Qm^o = 0.28 and h = 0.7 



Quantity 



ACDM 



EdS 



h 

dc 
Rp 
Mh 

Halo profile 



0.28 
0.7 

15 /i"^ Mpc 
1.1 h-^ Mpc 
1.5 • lO^^h-^MQ 
gaussian 
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FIG. 9: Shown are the magnification PDFs for the ACDM 
model described in Table Hill for a source at z = 1.5 in the 
case of No = 1 and No = 20 observations. 



and the EdS model with Qm,o = 1 and h = 0.5. The 
latter has been given a lower Hubble constant in order to 
better agree with the CMB observations. For simplicity, 
we are again taking our objects to have a constant proper 
radius Rp^ although for large, possibly unvirialized struc- 
tures this is not necessarily true. We shall come back to 
this issue in a forthcoming publication [62 . The size of 



the halos in Table HI roughly corresponds to the scale of 
large superclusters and has been chosen in order to have 
an average density contrast of 200 at virialization which 
we assume to happen at z = 1.6. For cluster scale ha- 
los the most appropriate profile to be used would be the 
NFW profile [58]. In this study, however, we have used 
for simplicity the gaussian profile which does not need 
any extra parameter. Again we will investigate various 
profiles in a forthcoming publication [62 . 

We plot the magnification PDF for the models of Ta- 



ble [Hi] in Figs. [9][T0] for the ACDM and for the EdS mod- 
els, respectively. As expected, the distributions are now 
significantly more skewed than they were in the case of 
the smaller halos of the previous Section: in both cases 
the fundamental PDF (No = 1) has the mode at the 
demagnification corresponding to the empty beam, the 
maximum demagnification possible. Moreover Figs. 9][TQ 
clearly illustrate the importance of introducing the dis- 
persion indicators of Eqs. (95). Even if the a relative to 
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FIG. 10: Same as in Fig. [9] but now for the EdS-model 
described in Table |III| Also displayed is the definition of 
FWHM- and FWHM+ used in Eqs. (|95). 



Am <C (Tm- If the model universe of Table [TTI| captures 
the actual degree of inhomogeneity of the universe, then 
these results suggest that lensing bias has to be incorpo- 
rated within the data analysis of SNe observations for the 
present-day data sets. Let us remind that the crucial fea- 
ture that our model universe captures, and which gives 
rise to the large biases, is that photons can travel through 
voids and miss the localised over densities. This feature is 
absent from swiss-cheese models where the boundaries of 
the holes are designed to have compensating overdensi- 
ties. In these models a photon that passes through a void 
always has to pass also through a compensating high- 
density shell, which results in a constrained PDF. It is 
not surprising then that such models have been shown to 
have on average little lensing effects [151 IS] • 



VI. CONCLUSIONS 



the PDF of No = 1 is larger than the a relative to the 
PDF of No = 20 because of the long high-magnification 
tail, the dispersion v behaves in the opposite way show- 
ing that with little data sets the actual dispersion is much 
smaller than the one computed assuming almost-gaussian 
distributions. 

In Fig. [TT] we show the magnification bias Am. Left 
panels correspond to the ACDM and right panels to the 
EdS model as the GBS. As we have explained, the bias of 
the effective PDF depends quantitatively on the number 
of observations at each redshift. The upper panels show 
the bias computed from No corresponding to the UC- 
data displayed in the upper panel of Fig. |4j The line has 
been made continuous by using a continuous intepolation 
for No{z) between the actual data points. The lower 
panels show the results for No corresponding to the the 
simulated data of the JDEM experiment. Also plotted 
are the dispersions V- and V-^ introduced in Eqs. (95): 



they give a quantitative estimation of the dispersion rel- 
ative to the plotted bias. We stress that V- and V-^ refer 
to the PDF relative to No measurements and therefore 
they are good indicators of the actual dispersion relative 
to the data set used. 

First note that the EdS model has larger lensing ef- 
fects (observe that in Fig. 11 the scales in the left and 
the right panels are not the same). This is clearly 
seen from the convergence for the empty beam given 
in Eq. (21): f<iE h'^ ^m,o and even if the EdS uni- 
verse has a lower /i, the larger value of 1^m,o is dominant: 

application of the latter amplification of lensing effects in 
the EdS model. 

Moreover, the Union Compilation of Ref. [3] does not 
have enough SNe measurements at high redshifts {z > 1) 
to suppress the lensing bias, while a JDEM-like survey of 
2000 SNe will be able to recover the GBS result within 
the intrinsic SNe dispersion of (Jm ~ 0.1 mag, that is. 



In this paper we have presented a new method to cal- 
culate the magnification probability distribution function 
(PDF) for a universe made of randomly distributed ha- 
los. The method is based on the weak lensing approx- 
imation and on generating stochastic configurations of 
halos along the line of sight, or along the photon geodesic 
from source to the observer. The basic physical feature 
incoporated by the method is the fact that underdensi- 
ties occupy most of the volume while most of the mass 
lies in overdense regions in the form of clusters and fil- 
aments. As a consequence the column density along a 
single geodesic is likely to be lower than the average den- 
sity of the global background solution. This is of poten- 
tial importance for the interpretation of the supernova 
observations, which are still probing much smaller angu- 
lar scales than the scale at which the homogeneity of the 
GBS is recovered. 

We derived a simple statistical formula by use of which 
one can easily compute the fundamental PDF for a single 
event and the effective PDF for quantities averaged over 
a number of observations No- Our formulae were ex- 
plicitly written for arbitrary mass distributions of spher- 
ical halos with arbitrary density profiles and they can be 
straightforwardly extended to include also different halo 
geometries (say thin cylinders to model filaments). We 
also showed how one can easily incorporate most selection 
effects into the formulae for the observable fundamental 
PDF. Along with this paper, we released the turboGL 
package [l], a simple and very fast Mathematica code to 
perform numerical simulations based on our model. The 
code will be continuously updated to incorporate more 
features (halo mass distributions, geometries, systematic 
biases etc) in the future. 

While our method can easily compute the PDF relative 
to arbitrary halo mass distributions, selection biases and 
halo profiles and evolutions, we focused for simplicity on 
the simplest configuration of one family of halos with no 
evolution and no selection biases in the numerical exam- 
ples shown in the present paper. More general results will 
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FIG. 11: Lensing bias as a function of redshift for the parameters of Table ^O] The left panels are for the ACDM model, while 
the right panels are for the EdS model. Moreover the top panels use the redshift distribution for SNe measurements given by 
the Union Compilation, while the bottom panels use the redshift distribution for a JDEM-like survey of 2000 SNe, see Fig. [2] 
Filled circles mark the values obtained from the numerical simulation and the dotted curves give the dispersions v- and , as 
indicated in the figure. 



be presented in a forthcoming paper [62^ . We checked the 
validity of our weak lensing approximation against exact 
results of light propagation in an inhomogeneous universe 
and the validity of our stochastic approach against the 
ray-tracing simulations of Ref. [12 (in a model universe 
with halos of mass M ^ lO^^/i^^M©). In both cases we 
were able to recover the main results within a few percent 
accuracy. 

In addition, we considered the biasing effect on the 
PDF due to very large structres, M ^ lO^^/i^^M©, 
whose existence is suggested by the large voids and fil- 
amentary structures seen both in the large-scale simula- 
tions and in the galaxy redshift surveys ^Ml EH EB. • We 
produced simulated PDFs for such universes both in the 
ACDM and in the Einstein-de Sitter (EdS) background 
model. We also produced the distributions for binned 
sets of observations and computed the biases and disper- 
sions for these effective PDF. Our results suggest that 
the lensing bias could affect the extraction of cosmologi- 
cal parameters from the current best data sets. However, 
we found that a JDEM-like survey should be able to re- 
move the lensing effects of even the largest imaginable 
structures assuming, of course, that no selection biases 
were present in the measured SNe. 



Indeed, in addition to computing the statistical bias, 
and perhaps beyond that, the most potential use of our 
method may be in computing the effect of different sets of 
selection biases on the observational magnification PDF. 
For example, the effects of rejecting outliers, or existence 
of zones of avoidance in the sky, possibly correlated with 
the densest structures, are easily incorporated. This issue 
will be matter of further investigation [62]. The method 
can be extended to include also strong lensing effects and 
eventually it will be interesting to incorporate also the 
redshift effects to our simulation package. 



APPENDIX A: PARTIALLY FILLED BEAMS IN 
EDS UNIVERSE 

Here we shall compare the analytic expressions for ex- 
act and cumulative weak lensing convergence functions 
for empty and partially filled beams in the EdS model. 
Given the redhsift dependence of the co-moving distance 

r(^) = ^(l-(l + ^)-V2) (Al) 
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we can compute the empty bin demagnification of a 
source at redshift z (r^ = r(z)) 



'^^(^^ = -2 — 1 



r) ao 



a{t{r)) 



1 + fl + z)"^/^ 
12-3 . : ln(l + z) . (A2) 



-1/2 



Note that ^e^z) does not depend on Hq. That is, it only 
depends on the redshift difference and not on the time 
of the observation. For a partiahy fihed beam, with the 
fining factor 



Phe 



the weak lensing result becomes simply 
KtJz) = {l-a)KE{z), 



(A3) 



(A4) 



as shown in Section |III A| Let us now compare this 
result to the exact formulae from Refs. [H |53l [63] . 
Zel'Dovich [53 was perhaps the first to recognize the 
importance of inhomogeneities on FLRW distances. He 
found the following exact result for the luminosity dis- 
tance along an empty beam (a line of sight emptied of 
any matter) in the EdS universe: 

DlAz) = ^ (1 + z? (l - (1 + z)-"^) ■ (A5) 

This result can be extended to the case of a partially 
filled beam [63j: 



kHo 



), (A6) 



where: 



k = (25 - 24a; 



,1/2 



(A7) 



with a given in Eq. ( |A3[ ). Note that for a = Eq. ( |A6| ) 

reduces to Eq. (A5) and for a = 1 it reduces to the EdS 



luminosity distance: Dl,i{z) = (1 + z)aor(z), where r{z) 
is given by Eq. ( Al ). Using Eq. ( Al ) and Eq. ( |A6| ) we can 
then compute the (de) magnification /i = (l^L,a/^L,i)^ 
and the convergence: 



1/2 



1 - (^+^)^'"'^^' ^ (A8) 



-1/2 



Results ( A4) and (|A8|) look very different, but they agree 



numerically to within a few per cent up to z ~ 2, as can 
be seen from the Fig. [2] For small z the agreement of 



(A4) and (A8) can be seen explicitly from their power 



series expansions: 



and 



^z^)^0{z') (A9) 



3(1 4 
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z^ + (9(z^). (AlO) 



The first two terms of the expansions vanish in accor- 
dance with the fact that lensing is negligible at small 
redshift s. Moreover, the functions tz^ and tv^^^ agree 
exactly up to third order and up to the fifth order the 
correction is ~ (1 — a)^. That is, while for a ^ 1 the 
convergence goes to zero as 1 — a, the difference between 
the exact result and the weak lensing approximation goes 
to zero even faster (1 — a)^. 



APPENDIX B: IDENTITY OF EQ. (36) 



We want to express the following integral / as a volume 
integral of cp: 



/ = 27r / dhhVih.t) = 27r / dhh 

Jo Jo Jb 



2ldl 



First we rewrite the latter as 
cR pR 



pU Pit 1 

I = 27T db — (-Vl^-b^) 2ldlip{l,t). (Bl) 
Jo Jh do \ J 



Then we introduce the auxiliary function 

pR 



Z{h) 



Vl^ -b^2ldlip{l,t), 



and calculate its derivative: 



Z\b) = y/W^2ldlip{l,t) 



l=b 



db 



-b^2ldlip{l,t) 



+ ^"^^ (-V^2_52^ 2ldl^{l,t). 



Finally we substitute the latter expression in Eq. (Bl): 

pR 

I = 2iT dbZ\b) = 27r{Z{R)-Z{0)) 
Jo 

pRit) 

= 47r / dbb^if{b,t). 
Jo 
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